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Abstract 

Background: The trend for the fabrication of electrical circuits with nanoscale dimensions has led to impressive progress in the 
field of molecular electronics in the last decade. However, a theoretical description of molecular contacts as the building blocks of 
future devices is challenging, as it has to combine the properties of Fermi liquids in the leads with charge and phonon degrees of 
freedom on the molecule. Outside of ab initio schemes for specific set-ups, generic models reveal the characteristics of transport 
processes. Particularly appealing are descriptions based on transfer rates successfully used in other contexts such as mesoscopic 
physics and intramolecular electron transfer. However, a detailed analysis of this scheme in comparison with numerically exact 
solutions is still elusive. 

Results: We show that a formulation in terms of transfer rates provides a quantitatively accurate description even in domains of 
parameter space where strictly it is expected to fail, e.g., at lower temperatures. Typically, intramolecular phonons are distributed 
according to a voltage driven steady state that can only roughly be captured by a thermal distribution with an effective elevated 
temperature (heating). An extension of a master equation for the charge-phonon complex, to effectively include the impact of off- 
diagonal elements of the reduced density matrix, provides very accurate solutions even for stronger electron-phonon coupling. 

Conclusion: Rate descriptions and master equations offer a versatile model to describe and understand charge transfer processes 
through molecular junctions. Such methods are computationally orders of magnitude less expensive than elaborate numerical simu- 
lations that, however, provide exact solutions as benchmarks. Adjustable parameters obtained, e.g., from ab initio calculations allow 
for the treatment of various realizations. Even though not as rigorously formulated as, e.g., nonequilibrium Green's function 
methods, they are conceptually simpler, more flexible for extensions, and from a practical point of view provide accurate results as 
long as strong quantum correlations do not modify the properties of the relevant subunits substantially. 
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Introduction 

Electrical devices on the nanoscale have received substantial 
interest in the last decade [1]. Impressive progress has been 
achieved in contacting single molecules or molecular aggre- 
gates with conducting or even superconducting metallic leads 
[2,3]- The objective is to exploit nonlinear transport properties 
of molecular junctions as the elementary units for a future mole- 
cular electronics. While the initial experiments were operated at 
room temperature, low temperatures down to the millikelvin 
range, the typical regime for devices in mesoscopic solid state 
physics, are also accessible (see, e.g., [4-6]). This allows for 
detailed studies of phenomena such as inelastic charge 
transfer due to molecular vibrations [7-9], voltage driven 
conformational changes of the molecular backbone [10], 
Kondo physics [11], and Andreev reflections [6], to name but a 
few. 

These developments have been accompanied by efforts to 
advance theoretical approaches in order to obtain an under- 
standing of general physical processes on the one hand and to 
arrive at a tool to quantitatively describe and predict experi- 
mental data. For this purpose, basically two strategies have been 
followed. One is based on ab initio schemes that have been 
successfully employed for isolated molecular structures, such 
as, e.g., density functional theory (DFT). Combining DFT with 
nonequilibrium Green's functions (NEGF) allows us to capture 
essential properties of junctions with specific molecular struc- 
tures and geometries [2,3,12,13]. This provides insight into the 
electronic formations of contacted molecules and gives at least 
qualitatively correct results for currents and differential conduc- 
tances. However, a quantitative description at the level of accu- 
racy known from conventional mesoscopic devices still seems 
to be out of reach. Furthermore, these methods are not able to 
capture phenomena resulting from strong correlations effects, 
such as Kondo resonances. 

Thus, an alternative route, mainly inspired by solid state 
methodologies, starts with simplified models that are assumed 
to cover the relevant physical features. The intention then is to 
reveal fundamental processes characteristic for molecular elec- 
tronics that give a qualitative description of observations from 
realistic samples, but provide also the basis for a proper design 
of molecular junctions to exploit these processes. Information 
about specific molecular set-ups appears merely in the form of 
parameters which offer a large degree of flexibility. In general, 
to attack the respective many body problems, perturbative 
schemes have been applied, the most powerful of which are 
nonequilibrium Green's functions [14,15]. However, conceptu- 
ally simpler, easier to implement, and often better at revealing 
the physics, are treatments in terms of master or rate equations. 
Being approximations to the NEGF frame in certain ranges of 



parameters space, they sometimes lack the strictness of pertur- 
bation series, but have been extensively employed for meso- 
scopic devices [16] and quantitatively often provide solutions of 
at least similar accuracy. Roughly speaking, these schemes 
apply as long as quantum correlations between relevant subunits 
of the full compound are sufficiently weak [15]. Physically, it 
places charge transfer through molecular contacts in the context 
of inelastic charge transfer through ultrasmall metallic contacts 
(dynamical Coulomb blockade [17]) and in the context of 
solvent or vibronic mediated intramolecular charge transfer 
(Marcus theory) [18-20]. 

While rate descriptions have been developed in a variety of 
formulations before [21-28], the performance of such a frame- 
work in comparison with numerically exact solutions has not 
yet been addressed. The reason for this is simple: A numerical 
method that provides numerically exact data in most ranges of 
parameters space (temperature, coupling strength, etc.) has only 
very recently been successfully implemented in the form of a 
diagrammatic Monte Carlo approach [29] . Path integral Monte 
Carlo methods have been used previously for intramolecular 
charge transfer in complex aggregates [18,19] in a variety of 
situations, including correlated [30] and externally driven 
transfer [31] and, of particular relevance to the present work, 
transfer in the presence of prominent phonon modes [32], 

The goal of the present work is to study a simple yet highly 
nontrivial set-up, namely, a molecular contact with a single 
molecular level coupled to a prominent vibronic mode (phonon) 
which itself may or may not be embedded in a bosonic heat 
bath. We develop rate descriptions of various complexity, place 
them into the context of NEGF, and compare them with exact 
solutions. The essence of this study is, astonishingly enough, 
that rate theory provides quantitatively accurate results for mean 
currents over a very broad range of parameter space, even in 
domains where they are not expected to be reliable. 

Results and Discussion 

In subsection 1 we define the model and the basic ingredients 
for a perturbative treatment. A formulation which closely 
follows the P(E) theory for dynamical Coulomb blockade is 
discussed in subsection 2. Nonequilibrium effects in the 
stationary phonon distribution are analyzed in subsection 3 
based on a dynamical formulation of charge and phonon 
degrees of freedom. The presence of a secondary bath is incor- 
porated in subsection 4 together with an improved treatment of 
the molecule-lead coupling, which is exact for vanishing 
electron-phonon interaction. The comparison with 
numerically exact data and a detailed discussion is given in 
subsection 5. 
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1 Model 

We start with the minimal model of a molecular contact 
consisting of a single electronic level (dot) coupled to fermi- 
onic reservoirs, where a prominent internal molecular phonon 
mode interacting with the excess charge is described by a 
harmonic degree of freedom (Figure 1) [15,33,34]. 




Figure 1 : Single charge transfer through a molecular contact 
consisting of a single electronic level coupled to a harmonic phonon 
mode and contacted to metallic leads. Forward (no prime) and back- 
ward (with prime) rates are the basic ingredients for the approximate 
treatment, see text for details. 



Neglecting spin degrees of freedom the total compound Hamil- 
tonian is thus described by 



(spin-boson model [20]). Depending on the charge state of the 
dot the phonon mode is subject to potentials V q (x 0 ) = (ni(o 0 / 
2)(xq + log) 2 . Now, the presence of the leads acts (for finite 
voltages) as an external driving force alternately charging (q = 
1) and discharging (q = 0) the dot, thus switching alternately 
between Vq and V\ for the phonon mode. The classical energy 
needed to reorganize the phonon is the so-called reorganization 
energy A = V\{0) - V 0 (0) = Mq I h<o 0 . Quantum mechanically, 
the phonon mode may also tunnel through the energy barrier 
located around xrj = -/rj/2 separating the minima of Fnj. 

It is convenient to work with dressed electronic states on the dot 
and thus to apply a polaron transformation generating the shift 
Zq in the oscillator coordinate associated with a charge transfer 
process, i.e., 



U = exp -i 



(2) 



with momentum operator Pq = i^jhmaiQ 1 2 (b^ —b^) where b Q ' 
and /?o are creation and annihilation operators of the phonon 
mode, respectively. We mention in passing that complementary 
to the situation here, the theory of dynamical Coulomb blockade 
in ultrasmall metallic contacts is based on a transformation 
which generates a shift in momentum (charge) rather than pos- 
ition [17]. Now, the electron-phonon interaction is completely 
absorbed in the tunnel part of the Hamiltonian, thus capturing 
the cooperative effect of charge tunneling onto the 
dot and photon excitation in the molecule, i.e., 



H = H 



L/R 



-H D + H Fh - 



H T with 



H = H L/R + #T +// D +// D,Ph +// Ph 



: Z £ k,a c k,a c k,a + Z ( r k,a c k,a d + kc ^ 
a=L,R;k a=L,R; k (1) 



+e d cftd + — + —0- (x 0 + l 0 d f d) 
2m 2 V / 



Here, the !7\ a denote tun nel coupl ings between dot level and 
reservoir a and / 0 = Mq-^2 / hm^m contains the coupling Mq 
between excess charge and phonon mode. An external voltage V 
across the contact is applied symmetrically around the Fermi 
level such that = &o(k) + n a , with the bare electronic disper- 
sion relation sn(£) and chemical potentials jiy^ = +eVI2, //r = 
—eVI2. Below, this model will be extended further to include the 
embedding of the prominent mode into a large reservoir of 
residual molecular and/or solvent degrees of freedom acting as 
a heat bath. Qualitatively, since the dot occupation cftd can only 
take the values q = 0 or 1, the sub-unit Hd + H^ p^ + Hp^ 
describes a two state system coupled to a harmonic mode 



H Th =h(o 0 \ blb 0 +- 



(i \ (3) 

a=L,R; k \ n J 



Single charge tunneling through the device can be formally and 
exactly captured under weak conditions (e.g., instantaneous 
equilibration in the leads during charge transfer) within the 
Meir-Wingreen formulation based on nonequilibrium Green's 
functions [14,15]. For the current- voltage characteristics one 
finds 



I(V): 



4e 



jde 



L^R 



£l+£r 
x\iG > (e)-iG < (e) 



eV_ 
2 



8 + - 



eV 



(4) 
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with energy dependent lead self-energies 
Xa(s) = 27rXk!7k,al 2 S( s - s^) and with the Fourier transforms of 
the time dependent Green's functions G < {t) = i(d' d(t)) and 
G > (t) = -i(d(t)d^}. Upon applying the polaron transformation 
(Equation 2), one has 



G<( t) = i Jli> l \> m m 



(5) 



where all expectation values are calculated with the full Hamil- 
tonian (Equation 3). Of course, for 7^ a — > 0, the Green's func- 
tions factorize as, e.g., G < (t) -> i(d^d(t)} D exp[J(0] with the 
phonon correlation 



This is the adiabatic regime. In this latter range the impact of 
the adiabaticity on the diabatic ground state wave functions is 
weak for mo < 1 when the distance of the diabatic surfaces is 
small compared to the widths of the ground states. For mo > 1 in 
both regimes electron transfer is accompanied by phonon 
tunneling through energy barriers separating the minima of 
adiabatic or diabatic surfaces. The dynamics of the total com- 
pound are then determined by voltage driven, collective 
tunneling processes. Master equation approaches to be investi- 
gated below, rely on the assumption that both sub-units, charge 
degree of freedom and phonon mode preserve their bare phys- 
ical properties even in the case of finite coupling m 0 . Hence, 
since the model (Equation 1) can be solved exactly in the limits 
ntQ = 0 and o = 0 and following the above discussion, we expect 
them to capture the essential physics quantitatively in the 
domain mo < 1 and for all ratios o. We note that recently the 
strong coupling limit including the current statistics has been 
addressed as well [35,36]. 



Jit) 



:Po l o tPoWo \ 



(6) 



Ph 



into expectation values with respect to the dot (D) and the 
phonon (Ph), respectively. Any finite tunnel coupling induces 
correlations that in analytical treatments can only be incorpo- 
rated perturbatively. There, the proper approximative scheme 
depends on the range of parameter space one considers. Gener- 
ally speaking, there are four relevant energy scales Xl/R> Mq, 
k#T, and hmo of the problem corresponding to three inde- 
pendent dimensionless parameters, e.g., 



m 0 = , 9 = co 0 /ip , a = 



h(On 



(7) 



In the following we are interested in the low temperature 
domain 6 > 1 where thermal broadening of phonon levels is 
small such that discrete steps appear in the I-V characteristics. 
Qualitatively, seen from the dynamics of the phonon mode, two 
regimes can be distinguished according to the adiabaticity para- 
meter X//zff>o = o: For o < 1 the phonon wave packet fulfills, on 
a given surface V 0 or V\, multiple oscillations before a charge 
transfer process occurs. The electron carries excess energy due 
to the finite voltage, and this energy may be absorbed by the 
phonon to promote reorganization to the new conformation (in 
the classical case the reorganization energy A). In the language 
of intramolecular charge transfer this scenario corresponds to 
the diabatic regime with well-defined surfaces F q . In the oppo- 
site regime a > 1 charge transfer is fast such that the phonon 
may undergo multiple switchings between the surfaces Vq i. 



2 Rate approach I 

The simplest perturbative approach considers the cooperative 
effect of electron tunneling and phonon excitation in terms of 
Fermi's golden rule for the tunneling part H T . For this purpose 
one derives transition rates for sequential transfer according to 
Figure 1 . A straightforward calculation for energy independent 
self-energies Xl/R (wide band limit) gives the forward rate onto 
the dot from the left lead 



r L (K,8 D ): 



h 



eV 



^o( £ - £ d). 



(8) 



where /p(e) is the Fermi distribution. Inelastic tunneling asso- 
ciated with energy emission/absorption of phonons is captured 
by the Fourier transform of the phonon-phonon correlation 
exp[J(?)] leading to 



P 0 (E) = e-Pa-Pe^^[ s _ /ic0o(/ _i ) ] 
k,l 



(9) 



with Pgy e ={m\ 1 2)[coth(6 / 2) + 1] denoting the mean values 
for single phonon absorption (a) and emission (e). The exponen- 
tials in the prefactor contain the dimensionless reorganization 
energy m\ = A/h(»o- Apparently, inelastic charge transfer 
includes the exchange of multiple phonon quanta according to a 
Poissonian distribution. Further, one has the detailed balance 
relation Po(s) = e~P £ Po(s). For vanishing phonon-electron 
coupling m 0 — > 0 only the elastic peak survives, thus -Po( e ) — * 
8(s). We note again the close analogy to the P(E) theory for 
dynamical Coulomb blockade [17]. Moreover, golden rule rates 
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for intramolecular electron transfer between donor and acceptor 
sites coupled to a single phonon mode are of the same form 
with the notable difference, of course, that in this case one has a 
discrete density of states for both sites [20,22]. The funda- 
mental assumption underlying the golden rule treatment is that 
equilibration of the phonon mode occurs much faster than 
charge transfer. In the last two cases this is typically guaranteed 
by the presence of a macroscopic heat bath (secondary bath) 
strongly coupled to the prominent phonon mode. Here, the 
fermionic reservoirs in the leads impose phonon relaxation due 
to charge transfer only. Thus, for finite voltage the steady state 
is always a nonequilibrium state that can only roughly be 
described by a thermal distribution of the bare phonon system 
(see below). One way to remedy this problem is to introduce a 
phonon-secondary bath interaction as well (see below in 
subsection 4). The remaining transition rates easily follow due 
to symmetry 



r R (F,B D ) = -*r L (K,-E D ), 

rk(^s D ) = ^r L (-F,s D ), 
ri(K,s D ) = r L (-K,-s D ). 



(10) 



Now, summing up forward and backward events, the dot popu- 
lation follows from 



d Pdo 
dt 



" r tot,0 /'dot + r d 



(11) 



with the total rate r tot o = + + T'l + T'r and the rate for 
transfer towards the dot Tq = Tl + Fr obtained according to 
Equation 8. Note that for vanishing electron-phonon coupling 
Mo = 0 one has /zr totjU (Mo = 0) = £l + Zr- The steady state 
distribution p iot — > p^ ot = r D /r tot 0 is approached with relax- 
ation rate r tot q. For a symmetric situation £l = Zr w i m £ D = 0 
one shows that p^ ot = 1/2 independent of the voltage, while 
asymmetric cases lead to voltage dependent stationary popula- 
tions. The steady state current is given by I(V) = (e/2)[(r^ - 
T'r)(1 -Pdot) " (T'l - r R )pdot] such that 



I(V). 



r L r R -r' L r' R 

r tot,0 



(12) 



A transparent expression is obtained for 8q = 0, namely, 

eV 



i{vy 



2 



/p|e + — ||W(13) 



Despite its deficiencies mentioned above, the golden rule treat- 
ment provides already a qualitative insight into the transport 
characteristics. Typical results are shown in Figure 2. 
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Figure 2: /-V-characteristics for symmetric coupling Y.L = Zr ar, d f° r 
varying electron-phonon coupling mo at inverse temperature 6 = 25 
(solid) and 6=10 (dashed). 



The I—V curves display the expected steps at 
eV = 2nfiaiQ,n e Z . Each time the voltage eVI2 exceeds multi- 
ples of /icon new transport channels open associated with the 
excitation of one additional phonon. For higher temperatures 
the steps are smeared out by thermal fluctuations. The range of 
validity of this description follows from the fact that a factor- 
izing assumption for the electron-phonon correlation and an 
instantaneous equilibration of the phonon mode after a charge 
transfer has been used, which means that a < 1 and mq < 1 . The 
latter constraint guarantees that conformational changes of the 
phonon distribution remain small. 

There are now three ways to go beyond this golden rule approx- 
imation. With respect to the phonon mode, one way is to explic- 
itly account for the nonequilibrium dynamics, another is to 
introduce a direct interaction with a secondary heat bath in 
order to induce sufficiently fast equilibration. With respect to 
the dot degree of freedom one can exploit the fact that for 
vanishing charge-phonon coupling the model can be solved 
exactly. 

3 Master equation for nonequilibrated 
phonons 

To derive an equation of motion for the combined dynamics of 
charge and phonon degrees of freedom, one starts from a Liou- 
ville-von Neumann equation for the full polaron transformed 
compound (Equation 3). Then, applying a Born-Markov type of 
approximation with respect to the tunnel coupling to the fermi- 
onic reservoirs, one arrives at a Redfield-type equation for the 
reduced density matrix of the dot-phonon system [15]. An add- 
itional rotating wave approximation (RWA) separates the 
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dynamics of diagonal (populations) and off-diagonal (coher- 
ences) elements of the reduced density. Denoting with P q n the 
probability to find q charges on the dot (here, for single charge 
transfer q = 0,1) and the phonon in its f?-th eigenstate, one has 



dt 



I L/n,k 
a=L,R; k 



(14) 



/p(c)^;-/p(-c)^ 



with v 0 = 1, vi = -1 and energies E^=h(o Q (k - «)+v q (|j. a +8 D ). 
The matrix elements of the phonon shift operator 
/n,k = (" I ex P('Po l O lh)\k) read 



./n,l 



-'»o 



( „„ ~^ n ~k\ 



max(n,k) 

n ' 

^/=min(n,k+l) j 



(n-k)\ 

X\Fi ^max(«, k) + 1, | n — k \ +l,-m^j, 



(15) 



where iF\ denotes a hypergeometric function. The underlying 
assumptions of this formulation require weak dot-lead coupling 
o < 1 and sufficiently elevated temperatures ad < 1 for a 
Markov approximation to be valid. Although we will see below 
when comparing low temperature results with numerically exact 
solutions that this seems to be only a weak constraint. 

The calculation of the steady state distribution 
Pq = lim^^ Pq (t ) reduces to a standard matrix inversion. One 
can show that for a symmetric system with sq = 0, Y.L = Zr one 
has Pq = if 1 . A typical example for the mean phonon number 
( n ) = ^ qn n Pq is depicted in Figure 3. The curve is well 
approximated by a/mn with a ~ 0.7. Apparently, <n> diverges 
for niQ — > 0 since then Pq and P^ approach constants inde- 
pendent of the phonon number. Upon closer inspection one 
finds that excitation is more likely than absorption, i.e.,f(n,n + 
1) >f (n,n - 1), for all 0 < n < No(niQ) where TVo(wio) increases 
with decreasing rtiQ. The opposite is true for n > No(mo) such 
that in a steady state, depending on the voltage, the tendency is 
to have higher excited phonon states occupied by smaller 
couplings rriQ. In particular, for strong coupling transitions n — > 
n + k,k > 0 are blocked at small n. 
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Figure 3: Mean phonon number in nonequilibrium for eV- 
versus the electron-phonon coupling mo. 


3ficjQ and 



reliable, as Figure 4 reveals. While it clearly shows the general 
tendency of a substantial heating of the phonon degree of 
freedom induced by the electron transfer, the profile of a 
thermal distribution strongly differs from the actual steady state 
distribution. 



0.20 r 
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Figure 4: Phonon number distribution in nonequilibrium for eV = 
5fiu)o, mo = 0.5 and kgTlhuio = 0.1 (histogram). The solid line depicts 
a fit to a Boltzmann distribution. See text for details. 



Nonequilibrated phonons leave their signatures also in the I—V 
curves as compared to equilibrated ones. The net current 
through the contact follows from the summing up of the transfer 
rates from/onto the dot according to Equation 14, hence, 



A convenient strategy to include nonequilibrium effects in the 
phonon distribution, sometimes used in the interpretation of 
experimental data, is the introduction of an effective tempera- 
ture 7 e ff. This way one could return to the simpler modeling of 
the previous section. However, the procedure to identify 
P 0 n + Tf « P p n = exp(-p eff hv 0 n) I [1 - exp(-(3 eff ta 0 )] is not 



I{V)- 



2h 



n,k 



^L/p(4' k L )- 2 R/p( £ nk R ) 



nil 
^0 



(16) 



421 



Beilstein J. Nanotechnol. 2011, 2, 416-426. 



Figure 5 shows that deviations are negligible for low voltages in 
the regime around the first resonant step (\eV/2\ < h(£>o), where 
at sufficiently low temperatures only the ground state partici- 
pates such that the steady state distribution coincides with the 
thermal one. For larger voltages deviations occur with the ten- 
dency that for smaller couplings mo the nonequilibrated current 
is always smaller than the equilibrated one (/ non < ^eqX while 
the opposite scenario (7 non > 7 eq ) is observed for larger m 0 . At 
sufficiently large voltages, one always has 7 non < 7 e q- This 
behavior results from the combination of two ingredients, 
namely, the phonon distributions P„ and the Franck-Condon 
overlaps |/n,kl 2 - To see this in detail, let us consider a fixed 
voltage. Then, on the one hand, for smaller niQ the steady state 
distribution is broad (cf. Figure 3), such that, due to normaliza- 
tion, less weight is carried by lower lying states compared to a 
thermal distribution at low temperatures; on the other hand, for 
niQ < 1 the overlaps \f n ^\ 2 favor contributions from low lying 
states in the current (16), which is thus smaller than 7 e q. For 
increasing electron-phonon coupling mo > 1, the overlaps \f n j<\ 2 
tend to include broader ranges of phonon states also covered by 
-Pq" 1 , compared to those of low temperature thermal states. A 
voltage dependence arises since with increasing voltage higher 
lying phonon states participate in the dynamics supporting the 
scenario for smaller couplings. Interestingly, as already noted in 
[23] the overlaps \f n J may vanish for certain combinations of 
n,m depending on mg due to interferences of phonon eigenfunc- 
tions localized on different diabatic surfaces K q , where q = 0,1. 



subsection 2 to a situation where the secondary bath is charac- 
terized by Gaussian fluctuations. Its corresponding modes can 
thus effectively be represented by a quasi-continuum of 
harmonic oscillators for which the phonon correlation function 
(Equation 6) can be calculated easily 



J(t)= L — 



co dm 1(a) 



coth 



ahfi 
2 



cos(cof) - 1] - i sin(co?) 



(17) 



Here the spectral density 1(a) now describes the combined 
distribution of the prominent mode and its secondary bath. It is 
thus proportional to the imaginary part of the dynamic suscepti- 
bility of a damped harmonic oscillator [20]. For a purely ohmic 
distribution of bath modes, one has 



2 3 

I ((a) = 2ffJ 0 © 0 



yea 



j 9 9 9 9' 

(or - cog) + y or 



(18) 



where y denotes the coupling between phonon mode and bath. 
The Fourier transform of exp(J) reads at finite temperatures 



P y (e) = 



; p T5(e) 



WW) 



k 



Py,e 



(19) 



k\ H s + hClk-hCll 
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Figure 5: /-(/-characteristics for equilibrated (solid) and nonequili- 
brated (dotted) phonon distributions according to Equation 13 and 
Equation 16, respectively. 



4 Rate approach II 

The assumption of a thermally distributed phonon degree of 
freedom during the transport can be physically justified only if 
this mode interacts directly and sufficiently strongly with an 
additional heat bath (secondary bath) realized, e.g., by residual 
molecular modes. Here we will generalize the formulation of 



with the frequency Q given by Q. = con^ + iy/2 and 
p y a (Q) = (ml 1 2^Q 2 )[coth(pftQ / 2) - 1] where the parameter % 
is £ j = ^j\-y 2 /4coo . Further, py e (£l) = -py e (Q*) (* means 
complex conjugation) and p T = 9?[py ifl + py, e ]/2. In the above 
expression, contributions from the Matsubara frequencies in 
Equation 17 have been neglected, since they are only relevant in 
the regime yh$ » 2n, which is not studied here. Apparently, 
the coupling to the bosonic bath effectively induces a broad- 
ening of the dot levels hy(k + 1)12 compared to the purely elastic 
case (Equation 9). In the low temperature regime, where for 
equilibrated phonons absorption (related to k) is negligible, the 
widths grow proportionally to /. The presence of the secondary 
bath drives the prominent phonon mode towards thermal equi- 
librium with a rate proportional to this broadening. Hence, if the 
time scale for thermal relaxation is sufficiently smaller than the 
time scale for charge transfer, i.e., 1/xi = (£l + ZrVy <<c 1> the 
assumption of an equilibrated phonon mode is justified and the 
golden rule formulation (Equation 13) can be used with 
Pn(s) — * Py(e). However, this argument no longer applies in the 
overdamped situation y/ojQ » 1, where the phonon mode 
exhibits a sluggish thermalization on the time scale y/coq, 
which may easily exceed x\. 
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As already mentioned above, for vanishing charge-phonon 
coupling niQ = 0, the model (Equation 3) can be solved exactly 
for all orders in the lead-dot coupling [15]. In the frame of a 
rate description, one observes that in this limit the dot popula- 
tion (Equation 11) decays proportionally to (£l + Zr). The 
golden rule version of the theory neglects this broadening in 
Equation 13 since it is associated with higher order contribu- 
tions to the current (Equation 13). Now, recalling that Po(s) 
reduces to a delta function for mo — > 0, this finite lifetime of the 
electronic dot level is included for all orders by performing the 
time integral in the Fourier transform with s — > s - ;'(Xl + ZrV 
2 = s- i'r tot (Mo = 0)/2 [see Equation 11]. In fact, this way one 
reproduces the exact solution (one electronic level coupled to 
leads with energy independent couplings), i.e., its exact spec- 
tral function. To be specific, let us restrict ourselves for the 
remainder of this discussion to the symmetric situation 
Zl = Zr = 272, and sd = 0. Then, in the presence of the phonon 
mode (m Q i= 0) the corresponding function P 0 £ (e) follows from 
Equation 9 by replacing the delta function by + h(o(k - I) + 
(£/2]. Again following the idea of a rate treatment, an impro- 
ved version of this result accounting for higher order 
electron-phonon correlations is obtained through the decay rate 
ftot.oO^o 0), instead of the bare dot level width J^/h = 
r tot; o(M 0 = 0). Equivalently, one replaces il[e + h(n(k - I) + C£J 
2] — > i'/[s + hm(k - I) + iT tot o/2] to arrive at an improved 
P^(s). We note that within a Green's function approach, and 
upon approximating the corresponding equations of motion, a 
similar result has been found in [15,33], with the difference 
though that instead of r tot; o an imaginary part of a phonon state- 
dependent self-energy Z£ | appears. One can show that the r tot 0 
appearing here within a rate scheme is related to a thermally 
averaged /j 2 . 

Now, an additional secondary bath can be introduced as above 
by combining Equation 19 with Pj Z , leading eventually to 



l(k,l)>(0,0) 



k 1 

Py,a Py,e 



(20) 



scheme is the following: It applies to all couplings a in the 
domain where the electron-phonon coupling is weak m 0 < 1 . In 
particular, second order processes in o capture cotunneling 
processes. For mo > 1 charge transfer is strongly suppressed and 
the phonon dynamics still occurs on diabatic surfaces for o « 1 
so that we expect the approach to cover this range as well. 

5 Comparison with numerically exact results 

A numerically exact treatment of the nonequilibrium dynamics 
of the model considered here is a formidable task. The number 
of formulations which allow simulations in nonperturbative 
ranges of parameter space is very limited. Among them is a 
recently developed diagrammatic Monte Carlo approach 
(diagMC) based on a numerical evaluation of the full Dyson 
series, which, in contrast to numerical renormalization group 
(NRG) methods [37], covers the full temperature range. For 
calculations of single charge transfer, results have been 
obtained with and without the presence of a secondary bath 
interacting with the dot phonon mode. 

We note that computationally these simulations are very 
demanding as for each parameter set and a given voltage the 
stationary current for the I-V curve needs to be extracted from 
the saturated value of the time dependent current I(t) for longer 
times. Typical simulation times are on the order of several days 
to weeks, depending on the parameter range. In contrast, rate 
treatments require minimal computational effort and can be 
done within minutes. Here, we compare numerically exact find- 
ings with those gained from the various types of rate/master 
equations discussed above. 

We start with the scenario where the coupling to a secondary 
bath is dropped (y = 0) to reveal the impact of nonequilibrium 
effects in the phonon mode. The formulation for an equili- 
brated phonon is based on Equation 13 with Pq replaced by P^j 
in Equation 20, while the steady state phonon distribution is 
obtained from the stationary solutions to Equation 14. In the 
latter approach the intrinsic broadening of the dot electronic 
level due to coupling to the lead is introduced in the following 
way: One first determines via Equation 14 a steady state distrib- 
ution P q . This result is used for an effective self-energy contri- 
bution (total decay rate) for nonequilibrated phonons, i.e., 



s + hQ.k -h€l 1 + ihT tot 0 1 2 



tot,neq 



(V)- 



Z l/n,k 
<x=L,R; n,k,q 



P q"/p (^nk)' 



(21) 



The width of the electronic dot level is thus voltage dependent 
and approaches the bare width from below for large voltages, 
that is limp_ >oo r tot; o(K) = YJh- The range of validity of this 



where = ha 0 (n - k) + n a . We note in passing that 
lim F _^ oo r totineq (F) = YJh = r to t(M 0 = 0). Subsequently, an im- 
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Figure 6: /-(/characteristics according to approximate models based 
on equilibrated phonons (solid) and nonequilibrated phonons (dashed) 
together with exact DQMC data (dots) for kgTIY. = 0.2 and without 
coupling to a secondary bath (y = 0). All quantities are scaled with 
respect to the dot-lead coupling £. 



proved result for the steady state phonon distribution at a given 
voltage is evaluated working again with Equation 14, but using 
the replacement: 



fj h(O 0 (k-l)T 



eV 



£=F- 



eV 



hr 



tot,neq 



(22) 



2 J[£-ho) 0 (k-l)] 2 +h 2 T 2 



tot,neq 



Of course, for £ — > 0 the standard Fermi distribution is 
regained. The corresponding steady state phonon distribution 
eventually provides the current according to Equation 16 with 
the same replacement (Equation 22) in this expression. The 
procedure relies on weak electron-phonon coupling mq < 1 and 
in principle also requires sufficiently elevated temperatures. 



Results are shown in Figure 6 together with corresponding 
diagMC data for various coupling strengths thq. Interestingly, 
the equilibrated model describes the exact data very accurately 
from weak up to moderate electron-phonon coupling mo ~ 1, 
while deviations appear for stronger couplings m 0 > 1 and volt- 
ages beyond the first plateau eV> 2 /icon. For win > 1 nonequilib- 
rium effects are stronger and the corresponding master equa- 
tion (Equation 14) gives a better description of higher order 
resonant steps. Moreover, as already addressed above, even in 
this low temperature domain the approximate description 
provides quantitatively reliable results. 

In Figure 7 the frequency of the phonon mode is fixed and only 
the electron-phonon coupling is tuned over a wider range. For 
strong coupling (here mq = 2) the equilibrated (nonequilibrated) 
model predicts a smaller (larger) current than the exact one in 
contrast to the situation for smaller thq. This phenomenon 
directly results from what has been said above in subsection 3: 
For stronger coupling the Franck-Condon overlaps favor higher 
lying phonon states that are suppressed by a thermal distribu- 
tion. 

After all, the approximate models give not only a qualitatively 
correct picture of the exact I—V curves, but even provide a rea- 
sonable quantitative description in this low temperature domain. 

In a next step the coupling to a secondary bath is turned on 
(y i= 0) enforcing equilibration of the phonon mode, see Equa- 
tion 20. The expectation is that in this case departures from the 
equilibrated model are reduced. In Figure 8 data are shown for a 
ratio niQ = 4/5, where deviations occur at larger voltages, as 
observed in the previous figures. Obviously, due to the damping 
of the phonon mode the resonant steps are smeared out with 
increasing y. However, the approximate model predicts this 
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Figure 7: As Figure 6 but for fixed fiuVX = 5 and varying 
electron-phonon coupling. 
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interacting with the phonon with various coupling constants y. Shown 
are approximate results (solid) using Equation 20 and diagMC data 
(dots); energies are scaled with X- Other parameters are kg,TIY = 0.2, 
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effect to be more pronounced as compared to the exact data, 
particularly for stronger coupling hyf£ > 1, while still y/coo < 1. 
In fact, in the limit of very large coupling only the k = I = 0 
contribution to Equation 20 survives, such that at zero tempera- 
ture one arrives at 



( 



lim I{V) = 4 



- arctan 



eV 



V 



(23) 



with the current at large voltages I x = eYJ^h and r^oC 7 ) < 
£//z, where equality is approached for V — > oo. It seems that a 
broadened equilibrium distribution of the phonon, induced by 
the secondary bath according to Equation 20, overestimates the 
broadening of individual levels. Since the approach is exact in 
the limit mo — > 0, the deviations appearing in Figure 8 are due 
to intimate electro n-phonon/secondary bath correlations not 
captured by the rate approach. In the overdamped regime, i.e., 
y/cog > 1, the dynamics of the phonon mode slow down and may 
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Figure 9: As Figure 6, but for nonequilibrated phonons based on an 
extended master equation (solid) in comparison to exact diagMC data 
(dots). 
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become almost static on the time scale of the charge transfer. In 
this adiabatic regime an extended version of the master equa- 
tion (Equation 14) is not trivial since the conventional eigen- 
state representation becomes meaningless. It would be better 
then to switch to phase-space coordinates and develop a formu- 
lation based on a Fokker-Planck or Smoluchowski equation for 
the phonon. This will be the subject of future research. 

The essence of this comparison is that, as anticipated from 
physical arguments already in subsection 1, a rate description 
does indeed provide quantitatively accurate results in the regime 
of weak to moderate electron-phonon coupling mo < 1 and for 
all o. Deviations that occur for larger values of mq can partially 
be explained by nonequilibrium distributions in the phonon 
distribution, where, however, the master equation approach 
seems to overestimate this effect. In order to obtain some 
insight into the nature of this deficiency, a minimal approach 
consists of extending Equation 14 with Equation 22, a mecha- 
nism that enforces relaxation to thermal equilibrium with a 
single rate constant Tq that serves as a fitting parameter. 
Accordingly, the respective time evolution equation for P„ (f) 
receives an additional term -ro[Pq(t) - Pp] with the Boltz- 
mann distribution for the bare phonon degree of freedom Pp" . 
Corresponding results for the same parameter range as in 
Figure 6 and Figure 7 are shown in Figure 9 and Figure 10 
including comparison with the exact diagMC data. There, the 
same equilibration rate HTqIY. = 0-25 is used for all parameter 
sets. Astonishingly, this procedure provides excellent agree- 
ment over the full voltage range. It improves results particu- 
larly in the range of moderate to stronger electron phonon 
coupling, but has only minor impact for niQ < 1. The indication 
is thus that electron-phonon correlations neglected in the orig- 
inal form of the master equation have effectively the tendency 
to support faster thermalization of the phonon. Indeed, prelimi- 
nary results with a generalized master equation, where the 
coupling between diagonal (populations) and off-diagonal 
(coherences) elements of the reduced charge-phonon density 
matrix is retained (no RWA approximation), indicate that this 
coupling leads to an enhanced phonon-lead interaction and thus 
to enhanced phonon equilibration. 
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